

set matsize 11000
clear all

cd \Data

********************************************************************************
** The following do file estimates temperature response functions by vintage
** using monthly billing data spanning 2008-2013. The estimates from these
** temperature response functions are displayed in Figures A8 and A7.



** The following subsection of the code creates the daily values for the
** temperature spline during 2012 and 2013 on days with average daily
** temperatures above 62 deg F. These matrices will be used to predict the
** cooling energy usage that would occur during 2012-2013 using models fitted
** to the 2008-10 and 2011-13 billing data.

local thresh1=52
local thresh2=62
local thresh3=72

use cons_data.dta, clear
collapse (mean) temp_avg, by(year month day)
keep temp_avg
mkspline temp_avg_1 `thresh1' temp_avg_2 `thresh2' temp_avg_3 `thresh3' temp_avg_4 = temp_avg
keep if temp_avg_3 > 0
keep temp_avg_3 temp_avg_4
gen ones = 1

mkmat temp_avg_3 temp_avg_4, matrix(Temp)
mkmat ones, matrix(ones)


** The following subsection of the code simultaneoulsy estimates two temperature response functions
** for each vintage of homes. The first temperature response function uses billing data
** from period 1 (2008-10). The second temperature response function uses billing data
** from the same premises during period 2 (2011-13).

forvalues y = 1960/2004 {

	use billing_dataset.dta, clear
	sort premise_id
	merge m:1 premise_id using premise_characteristics.dta
	keep if year_built == `y'
	mkspline temp_1 `thresh1' temp_2 `thresh2' temp_3 `thresh3' temp_4 = temp_avg
	
	xtreg kwh period#(c.temp_1 c.temp_2 c.temp_3 c.temp_4) i.period, fe cluster(premise_id)
	matrix define beta_`y' = e(b)'
	matrix define var = e(V)
	matrix define var_`y' = vecdiag(var)'
	
	}

matrix define beta = (beta_1960)
matrix define var = (var_1960)
forvalues y = 1961/2004 {
	matrix define beta = (beta, beta_`y')
	matrix define var = (var, var_`y')
}

	
forvalues y = 1960/2004 {
	
	matrix define coef_`y'_per1 = (beta_`y'[5,1], beta_`y'[7,1])  /* These row vectors report the average response among homes built in year `y' to a one degree increase in temperatures on days with average temperature between 62-72 and above 72 during the 2009-10 bills */
	matrix define coef_`y'_per2 = (beta_`y'[6,1], beta_`y'[8,1])  /* These row vectors report the average response among homes built in year `y' to a one degree increase in temperatures on days with average temperature between 62-72 and above 72 during the 2011-13 bills */
	
	matrix define var_`y'_per1 = (var_`y'[5,1], var_`y'[7,1])  /* These row vectors report the variance of the average response among homes built in year `y' to a one degree increase in temperatures on days with average temperature between 62-72 and above 72 during the 2008-10 bills */
	matrix define var_`y'_per2 = (var_`y'[6,1], var_`y'[8,1])  /* These row vectors report the variance of the average response among homes built in year `y' to a one degree increase in temperatures on days with average temperature between 62-72 and above 72 during the 2011-13 bills */
	
}

matrix define coef = (coef_1960_per1, coef_1960_per2)
matrix define var = (var_1960_per1, var_1960_per2)
forvalues y = 1961/2004 {
	matrix define coef = (coef \ coef_`y'_per1, coef_`y'_per2)
	matrix define var = (var \ var_`y'_per1, var_`y'_per2)
}

forvalues y = 1960/2004 {
	forvalues p = 1/2 {
		matrix define cooling_`y'_per`p' = Temp * coef_`y'_per`p'
		matrix define cooling_`y'_per`p' = ones' * cooling_`y'_per`p'
	}

	matrix define cooling_`y' = (cooling_`y'_per1, cooling_`y'_per2)
}	

matrix define cooling = (cooling_1960)
forvalues y = 1961/2004 {
	matrix define cooling = (cooling \ cooling_`y') /* Each row displays the predicted aggregate 2012-13 cooling that would occur in a home built during year `y'. The first column displays predictions stemming form the 2008-10 temp response functions and the second column displays the predictions stemming from the 2011-13 temp response functions.*/
}


clear
set obs 45
gen year = _n

svmat coef
rename coef1 coef_per1_med
rename coef2 coef_per1_high
rename coef3 coef_per2_med
rename coef4 coef_per2_high

svmat var
rename var1 var_per1_med
rename var2 var_per1_high
rename var3 var_per2_med
rename var4 var_per2_high

svmat cooling
rename cooling1 cooling_per1
rename cooling2 cooling_per2


********************************************************************************
**Figure A8: Predicted 2012-13 Cooling Using 2008-10 & 2011-13 Billing Data

** Left Panel
scatter cooling_per1 cooling_per2 year

** Right Panel
gen diff = cooling_per2 - cooling_per1
scatter diff year


********************************************************************************
**Figure A7: Coefficients from 2008-10 & 2011-13 Temp Response Functions

** Top Panel
gen std_per1_med = sqrt(var_per1_med)
gen std_per2_med = sqrt(var_per2_med)
scatter coef_per1_med coef_per2_med year

** Bottom Panel
gen std_per1_high = sqrt(var_per1_high)
gen std_per2_high = sqrt(var_per2_high)
scatter coef_per1_high coef_per2_high year

